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A stochastic model for the dynamics of enzymatic catalysis in explicit, effective solvents under phys- 
iological conditions is presented. Analytically-computed first passage time densities of a diffusing 
particle in a spherical shell with absorbing boundaries are combined with densities obtained from 
explicit simulation to obtain the overall probability density for the total reaction cycle time of the 
enzymatic system. The method is used to investigate the catalytic transfer of a phosphoryl group in 
a phosphoglycerate kinase- ADP-bis phosphoglycerate system, one of the steps of glycolysis. The di- 
rect simulation of the enzyme-substrate binding and reaction is carried out using an elastic network 
model for the protein, and the solvent motions are described by multiparticle collision dynam- 
ics, which incorporates hydrodynamic flow effects. Systems where solvent-enzyme coupling occurs 
through explicit intermolecular interactions, as well as systems where this coupling is taken into 
account by including the protein and substrate in the multiparticle collision step, are investigated 
and compared with simulations where hydrodynamic coupling is absent, ft is demonstrated that the 
flow of solvent particles around the enzyme facilitates the large-scale hinge motion of the enzyme 
with bound substrates, and has a significant impact on the shape of the probability densities and 
average time scales of substrate binding for substrates near the enzyme, the closure of the enzyme 
after binding, and the overall time of completion of the cycle. 



^ I. INTRODUCTION 

PQ Biochemical reactions in the cell are often carried out 

Q through complex chemical networks consisting of many 
• ^ coupled elementary component steps"'^ . Even the eluci- 
""^ dation of the molecular-level mechanism which under- 
^^ lies the operation of a single component in such net- 
'~~' works is often a difficult task. Computer simulation is 
. playing an increasingly important role in such mecha- 

^ nistic studies but direct simulation of many biochemical 
f^ processes is challenging because they occur on a diverse 
C^ range of scales. This fact has prompted the development 
r of coarse-grain or mesoscopic methods that allow one to 

^^ circumvent some of the difficulties related to dynamics 
fT^ that take place on long space and time scales^ ^ . In en- 
^^ zyme kinetics long times scales can arise from the dif- 
Cn fusive approach of the substrate to the enzyme and the 

\ \ conformational changes in the enzyme in the course of 
^ the catalytic reactions it carries out. There have been nu- 
k> merous simulation studies of the effects of diffusion on en- 

rS zyme kinetics.^ ^. In this paper we describe how one may 

j^ construct a mesoscopic model of an enzymatic cycle that 
incorporates the diffusive approach of substrates to the 
enzyme based on the solution of the diffusion equation, 
along with a particle-based description of the enzymatic 
reaction that involves protein conformational changes, re- 
lease of the product and the return of the protein to its 
original conformation. 

The method is used to investigate a specific enzymatic 
reaction, 



ysis network. In particular, we focus on the forward reac- 
tion that involves the transfer of a phosphoryl group from 
1,3-bisphosphoglycerate (bPG) to ADP by the PGK en- 
zyme to form 3-phosphoglycerate (PG) and ATP. (Often 
it is the reverse reaction that is studied experimentally 
due to the instability of bPG;-) Phosphoglycerate kinase 
is a monomeric protein of moderate size (416 amino acid 
residues in the human isozyme studied here) found in 
all living organisms, with a highly conserved amino acid 
sequence across different life forms. Its structure, con- 
sisting of two equal-sized domains labeled by the N- and 
C-termini of the protein, is well-adapted to selectively 
bind two substrates: bPG binds to the N-terminal, while 
the nucleotide substrates, MgATP or MgADP, bind to 
the C-terminal domain of the enzyme. Structurally, the 
N- and C-domains consist of a 6-stranded parallel beta- 
sheet surrounded by alpha helices (see Fig. IT]). 
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catalyzed by the enzyme phosphoglycerate kinase 
(PGK). This reaction is an important step in the glycol- 



FIG. 1: The open conformation of phosphoglycerate kinase 
showing the N- and C-terminal domains of the protein. 

The mechanism for the enzymatic reaction, which in- 



volves large hinge bending motions of the domains of the 
proteiipSH^, has been the subject of many kinetic stud- 
ig j9 | i3| ti5] rpjjg activ ity of th e enzyme requires both sub- 
strates to be boundfii'Mill When both substrates bind, 
the enzyme undergoes a large-scale hinge-bending con- 
formational change that brings the substrates close to 
one another to catalyze the dephosphorylation of bPG. 
In this "closed" conformation, the transition state is sta- 
bilized, lowering the free energy barrier for the transfer 
of a phosphoryl group. Upon transfer, the enzyme is 
forced into an open configuration and the PG and ATP 
products are released. 



We shall be concerned with the enzymatic activity of 
PGK under physiological conditions in the cell where the 
binding process is diffusion limited.'^ The binding pro- 
cess is well suited to be modeled as a two-step process 
in which first the substrates diffuse freely into a region 
near the enzyme, and then are drawn into the binding 
sites on the enzyme. Thus, it is reasonable to utilize 
a hybrid, stochastic procedure that combines analytical 
calculations with explicit simulation. The first step in 
the process of computing the distribution of time scales 
of the catalytic activity of the enzyme can be estimated 
by calculating first-passage times for the substrates mov- 
ing into the vicinity of the enzyme, while the second step 
requires a more detailed dynamical simulation due to the 
influence of the enzyme on the dynamics of the substrate. 
There have been simulations of the domain motions of 
PGK using a variety of methods .I^^I^S Given the large size 
of the protein and the long time scales of the motions, 
a full molecular dynamics simulation of the second step, 
which involves binding of the substrates to the enzyme 
in solution, the hinge-bending motion of the enzyme- 
substrate complex, followed by the reaction of the sub- 
strates and final release of products coupled with the re- 
opening of the enzyme, is computationally demanding. 
Gonsequently, we develop a coarse grain description of 
this part of the enzymatic cycle that is particle-based, 
includes enzyme, substrates and solvent molecules explic- 
itly and retains many features of full molecular dynamics. 



The outline of the paper is as follows. The two steps 
of the enzymatic reaction dynamics, diffusive approach 
of enzyme and substrate and substrate binding and reac- 
tion, are described first. The mesoscopic model for the 
protein, substrates and solvent, along with a description 
of the interaction potentials that control the binding of 
the bPG substrate to the active site and conformational 
changes in the PGK protein, are the topics of Sec. [ll] 
Section IIIII discusses the various time scales involved in 
the diffusive encounter between the substrates and the 
enzyme and shows how the relevant first-passage times 
can be computed analytically. The results of simulations 
of the dynamics are reported in Sec. |IV| while the conclu- 
sions of the study are summarized in Sec. IV] 



II. PROTEIN AND ITS CATALYTIC ACTIVITY 
IN SOLUTION 

We consider a system containing PGK enzymes, along 
with substrate and solvent molecules. The enzyme exists 
in open and closed forms and binding of both substrates 
is necessary for large-scale conformational changes to oc- 
cur.^^ 15 17 -^g suppose that the ADP substrate is bound 
to the enzyme and construct a coarse-grain model of the 
protein interacting with the bPG substrate in the pres- 
ence of solvent. As discussed below, under physiological 
conditions, ADP binds quickly and the rate of the en- 
zymatic reaction is determined by the binding of bPG. 
The model of the enzymatic activity of PGK entails a 
description of the interactions of bPG with the enzyme 
as it binds to the active site, the conformational changes 
in the protein that lead to the reactive event and the re- 
lease of product and return of the protein to its original 
conformation. 



Netw^ork model of PGK and interactions w^ith 
substrate 



A coarse-grain network model of the PGK protein is 
constructed by replacing each amino acid residue with a 
single monomer bead and connecting the beads by links 
or bonds;^ 24-26 r^j^^ bound ADP substrate is treated as 
one of the protein beads, while the bPG substrate is also 
described in a coarse-grained fashion as a single bead. 
The set of bead coordinates R^'' — (Ri, R2, . . . , Rnp) 
specifies the configuration of the protein (P) and we let 
R denote the coordinate of bPG, henceforth called the 
substrate (S). The construction of the potential energy 
function that is responsible for the protein conforma- 
tional state and interactions between the protein and 
substrate are described in detail in Appendix A. Here we 
simply sketch the main elements that enter in the design 
of the potential function, yp5(R^'', R;<^), that is able 
to describe both conformational states of the protein, 
the binding of bPG to the active site, and the result- 
ing changes of protein conformational states that occur 
on substrate binding and product release^^. 

To construct a network model for PGK, protein 
database configurations built from crystallographic data 
were analyzed to determine a set of pairwise interactions 
between residues. Each of the 416 residues was repre- 
sented by a single monomer bead in a linear polymer 
representation of the protein, with the position of each 
bead taken to be the Cartesian coordinates of the alpha- 
carbon of the peptide. Both open and closed forms of 
the PGK molecule were taken from the initial and fi- 
nal protein database configurations generated from the 
morphing analysis of the confor matio nal change between 
open and closed conformationa^SESl jn the Database of 
Macromolccular Movementa^Sl. 

Pairs of beads separated by a distance r < 10 A were 
recorded, generating separate lists of indices for open and 



closed conformations. The interaction lists for open and 
closed configurations were then compared, and a set Be 
of common interaction pairs or links were identified and 
assigned bond potentials in the following way. For links 
in Be the bond length as well as the magnitude of the 
difference rco between the bond lengths in the open and 
closed conformations were computed. The links in Be 
were then grouped into two new subsets, Bhc and Bsc, 
containing hard (he) or soft common (sc) links, respec- 
tively, based on the value of the separation distance Vco, 
where links with r^o < 4 were identified as hard links. 
The list of common links was then compared with the 
lists of open links and closed links. This process yielded 
2891 hard-common links and 519 soft-common links. In 
this study, the ADP substrate is treated as a single bead 
that forms hard links with three different beads in the 
enzyme. 

Pairs that exist in either the list of open links or the list 
of closed kinks but not in both were sorted into soft-open, 
Bso, and soft-closed, B^x sets, respectively. There are 448 
soft-open links (so), and 619 soft-closed links (sx). 

Before the enzymatic reaction can occur, bPG must 
bind to the active site of the enzyme. The binding pocket 
of the enzyme for this substrate was defined by beads 
with coordinates (Rg,Ri,Rf), where Rq, Ri and R2 
are the coordinates of the alpha-carbon of the glycine 
residues 386, 387 and 388 in the amino acid sequence 
of the PGK enzyme. The binding interaction between 
the bPG substrate at position R and the enzyme was 
assumed to depend on both the distance between the 
substrate and the bead in the active site with coordi- 
nate RJ, \R — Ri\ — i?si, as well as the orientation 
of the substrate with respect to a coordinate frame de- 
termined by three beads defining the binding pocket of 
the enzyme. As the substrate binds it triggers conforma- 
tional changes in the protein that lead to hinge closing 
to bring the bPG and ADP substrates into proximity for 
the phosphoryl group transfer. Consequently as bPG in- 
teracts with the protein in the course of binding to the 
active site, the open protein configuration is destabilized 
with respect to the closed configurations, driving the en- 
zyme towards the closed conformation. To achieve this 
conformational change, the interaction potentials for the 
soft, non-common set of links were taken to depend on 
a reaction coordinate £,{Rsi), which is a function of the 
distance between bPG and the active site. The net effect 
of the combination of these contributions is a protein- 
substrate interaction potential, yps(R^^, R; ^), which 
can draw in the bPG substrate, bind it to the active 
site of the enzyme in the open configuration, and then 
cause the enzyme to undergo a conformational change 
from an open to closed configuration. The network model 
of the protein and the binding of the substrate to the 
open conformation leading to hinge closing is shown in 
Fig. [2] After binding has taken place, the phosphoryl 
group transfer reaction is carried out by treating the reac- 
tion coordinate ^ as an external control parameter whose 
value is determined probabilistically. When the reaction 




FIG. 2: (left) Open conformation of the network model of 
PGK showing the approach of bPG to the binding pocket 
of the enzyme, (right) Protein conformation after substrate 
binding has resulted in hinge closing to form the closed con- 
formation. 



is complete the closed configuration is unstable and the 
enzyme reopens, completing the cycle. 



B. Solvent and its interactions with the protein 
and substrate 



The system also contains Ng solvent molecules with 
positions, r^' = (ri, r2, . . . ,rArJ and velocities, v^'' = 
{vi,V2, ■ ■ ■ ,vnJ. The solvent evolution is modeled by 
multiparticle collision (MPC) dynamics.ISS In MPC dy- 
namics there are no intermolecular potentials among sol- 
vent molecules. Instead, solvent molecules propagate in 
the absence of solvent-solvent interactions and undergo 
multiparticle collisions at discrete times r that account 
for the effects of many real collisions during this time in- 
terval. More specifically, after the streaming step, solvent 
particles are assigned to cells with length £ for the pur- 
poses of carrying out multiparticle collisions. The center- 
of-mass velocity Vc of particles in a cell is computed for 
each cell c, and the velocities of the solvent particles rel- 
ative to the center-of-mass velocity are rotated around a 
randomly chosen axis by an angle chosen from a set of 
possible rotations. This "collision" step conserves linear 
momentum, energy and particle number, and is consis- 
tent with hydrodynamic flow^S^. The collision step for 
a particle i in cell c is therefore: 



+ a; • (vi -Vc), 



(2) 



where v'^ is the post-collision velocity of particle i and a; 
is a rotation matrix. 

When the system contains proteins and substrates dis- 
solved in the solvent, the evolution is described by hybrid 
molecular dynamics-multiparticle collision (MD-MPC) 
dynamics.'^ In such hybrid dynamics, while the solvent 
molecules interact among themselves through multipar- 
ticle collisions, they interact with the solute molecules 
through solvent-bead intermolecular forces, Vst- The to- 
tal potential energy of the system is therefore given by 
Vt = Vps + Vsb and Newton's equations of motion are 
used to evolve the system under this potential energy 
for time intervals t between MPC events. This hybrid 



dynamics also satisfies the conservation laws and cor- 
rectly describes hydrodynamic interactions among solute 
species and fluid flows in the solvent. 



Penetrating solvent model 

Hydrodynamic interactions and solvent dissipation can 
also be included in a heuristic way by dropping the direct 
interactions between solvent and solute molecules and 
instead including the solute beads into the MPC step 
of the dynamicspS In this scheme, the solvent particles 
evolve freely between collision steps while the coordinates 
and momenta of the enzyme and substrate are evolved 
through Newton's equations of motion under the Vps 
potential function. 

More specifically, to allow for interaction between the 
solvent and beads, the collision rule is modified to include 
the velocity of the beads in the local ccntcr-of-mass ve- 
locity of particles in a cell. The center-of-mass velocity 
is computed for a cell c containing N^ solvent particles of 
equal mass m and a single bead of mass M and velocity 
V via: 



M 
Mt 



NcTn 
Mt 



^«z, 



(3) 



where Mr = Nciri + M is the total mass of particles in 
the cell and Vi is the velocity of solvent particle i in cell 
c. The collision rule for the penetrating solvent model 
with hydrodynamics is defined as 






(4) 



for the bead velocity V and the solvent velocities Vi. 
Since the magnitude and direction are conserved in the 
rotation, particle number, linear momentum and energy 
are globally conserved, resulting in proper hydrodynamic 
flow. 



Penetrating solvent without hydrodynamics 

For the purpose of assessing the importance of hydro- 
dynamic interactions, it is useful to construct an alter- 
native model in which the hydrodynamic effects are not 
present .'^SElHaS The collision rule for the penetrating sol- 
vent model can be modified by defining the center-of- 
mass velocity of particles in a cell to be 



Mt 



Nsm 

1 

Mt 



(5) 



where Ns is drawn from a Poisson distribution with mean 
value pVc, where p is the number density of solvent in 
the system and Vc is the cell volume. The total mass is 
Mt = M + Ngin, and Vs is an effective solvent veloc- 
ity drawn from a Maxwell-Boltzmann distribution with 



mass Nsm. Since this velocity is drawn at each colli- 
sion step, the velocity of the solvent is uncorrelated from 
one collision step to another. In this model, explicit sol- 
vent particle dynamics is replaced by the action of the 
collision operator. Since the velocity of the fluid is com- 
pletely decorrelated after a single collision step, any dy- 
namic correlations associated with a small value of the 
ratio of the mean free path to cell length strictly vanish. 



III. ENZYMATIC CYCLE DYNAMICS 

Complete enzymatic cycles can be simulated using the 
mesoscopic dynamical scheme described in the previ- 
ous section. When the protein, substrates and solvent 
molecules are modeled as structureless particles, full MD- 
MPC dynamics has been used to study the effects of dif- 
fusion on enzyme kinetics.'^ However, in the conditions 
that pertain to the interior of a cell, even this multi- 
scale method will not be computationally efficient if both 
the internal dynamics of the enzymes and the diffusive 
motion are considered. Under physiological conditions 
the concentrations of both substrates in the cytoplasm 
are relatively small*" ^i (0.14 mM for ADP and 0.001 
mM for bPG), while the enzyme concentration is roughly 
0.1 mM. If the substrates and enzyme are uniformly dis- 
tributed in the volume, the radius of the spherical volume 
around the enzyme containing a single substrate molecule 
is roughly tadp = 142 A for ADP and rbpc = 734 A for 
bPG. The sphere containing a single enzyme has a ra- 
dius of rpGK ~ 158 A. Estimating the viscosity of the 
cytoplasm to be roughly 5 times that of water, namely 
77 = 0.005 Kg/(m-s), and assuming the substrates have 
an effective radius Rs ~ 5 A, the Stokes-Einstein law 

D = kBT/{6nr]Rs) gives a value of D = 910 A^/iis. 
Given these conditions, we shall see that the ADP sub- 
strate binds typically before 5 /is, whereas the binding 
time of the bPG is very broadly distributed over many 
decades and is the main factor determining the reaction 
time. For this reason we suppose that ADP is bound to 
the enzyme and focus on the binding of bPG. 

From these considerations it is evident that the en- 
zymatic dynamics has a significant diffusion-influenced 
component; therefore, it is computationally inefficient to 
follow individual trajectories of the diffusive dynamics of 
substrates and enzymes in the solvent for the long times 
needed for enzyme-substrate encounters. Consequently, 
it is useful to decompose the process into portions where 
the substrates diffuse in the solvent without directly in- 
teracting with proteins, and portions where these species 
interact through direct intermolecular forces. The diffu- 
sive portions of the dynamics can be treated to a good 
approximation by analytical methods, while in the inter- 
acting portions the mesoscopic dynamical scheme can be 
used to describe details of the binding, conformational 
changes and reaction. These considerations suggest a 
stochastic model for the cycle dynamics that combines 
these types of dynamical evolution. 



A. Stochastic model for enzyme dynamics 

Initially, suppose the bPG substrate moves diffusively 
in a volume with radius rbpc surrounding the enzyme 
without any influence on its motion due to the pres- 
ence of an enzyme. Since the concentration of enzyme 
is a factor of 100 times that of the bPG, the number 
of enzymes in this volume should be Poisson distributed 
with an average number of 100 enzymes in the volume 
if there is no correlation in the density of enzymes. We 
assume the binding of the bPG to any enzyme in the 
volume occurs in the following way: At any given time, 
the bPG is within the spherical volume with radius r2 of 
some enzyme, which is smaller than the volume around 
the enzyme that contains a single substrate molecule (see 
Fig. |3| . The substrate can either diffuse to the binding 
region of this enzyme, or out of its volume. The bind- 
ing probability is dependent on how far the substrate is 
from the enzyme. If the substrate diffuses out of the 
volume of the enzyme, the first passage time out of the 
spherical volume can be recorded. Subsequently, the po- 
sition of the bPG relative to another enzyme is assumed 
to be randomly distributed in the volume of this other 
enzyme, and the process is repeated until the substrate 
passes through the inner spherical volume of radius ri 
around an enzyme. 

The point where the substrate passes through the in- 
ner sphere is uniformly distributed on the surface of the 
sphere. After passing through the inner sphere, the sub- 
strate will either bind to the enzyme or move out of the 
inner sphere and pass through a sphere of intermediate 
size (with radius r, with ri < r^ <^ r2). Since the dy- 
namics of the substrate is influenced by the presence of 
the enzyme and the solvent flow around it, the dynamics 
is no longer diffusive and must be simulated explicitly as 
described in the previous section. Starting from a uni- 
formly chosen point on the surface of the sphere with 
radius ri, if the substrate does not bind to the active 
site, the particle continues to diffuse starting from a ra- 
dial distance of r^ and either will be reabsorbed by the 
inner sphere or pass out of the volume through the outer 
sphere. 

For most of the dynamical evolution, the substrate dif- 
fuses freely without explicit solvent flow effects or influ- 
ence from the enzyme. For this type of dynamics, analyt- 
ical solutions to the diffusion equation can be used. The 
final regime to be described consists of the dynamics of 
the substrate from the surface of the inner sphere with 
radius ri to the active site on the enzyme in the presence 
of solvent. This final regime should be simulated directly, 
since the hydrodynamic motion of the solvent influences 
both substrate and enzyme motion. 

More speciflcally, the algorithm can be stated as fol- 
lows: 

(1) At the initial time if the substrate is at r2 a position 
r, ri < r < r2, is randomly selected. 

(2) Given a uniformly distributed random number 
£,r € [0, 1], if ^r < -Pi(^) the substrate is absorbed at the 




FIG. 3; Structure of the model. For the system considered 
here, we have chosen r2 = 31.6, ri = 9, and ri = 7 in simu- 
lation cell length units. This choice of radial distances allows 
one to minimize the amount of numerical simulation required 
while allowing for good statistics for various numerically com- 
puted densities. 



ri boundary, otherwise it is absorbed by the r2 bound- 
ary. Here Pi{r) is the probability that the substrate is 
absorbed at the ri boundary in the infinite time limit. 

(3) If it is absorbed at r2, a time is drawn from 
P2{t\r), the first-passage time density for absorption onto 
a sphere with radius r2 starting a distance r from center, 
and used to update the cycle time. 

(4) If it is absorbed at ri, a time is drawn from 
Pi(t|r), the first-passage time density for absorption onto 
a sphere with radius ri starting a distance r from cen- 
ter, and used to update the cycle time. Starting at ri, 
a full mesoscopic dynamical simulation is then carried 
out until reaction occurs or the substrate reaches the r^ 
boundary. If the dynamics results in a reaction, the time 
for this to occur is added to the cycle time and the enzy- 
matic cycle is complete. If instead the substrate reaches 
ri without reaction, this time is added to the cycle and 
we return to step (2) to continue the dynamics until the 
cycle is complete. The boundary at r^ is chosen to be 
significantly larger than ri to minimize the blocking ef- 
fect of the enzyme leading to a non-uniform distribution 
of points of absorption on the absorbing sphere. The ex- 
plicit forms of the Pi,2(?') and Pi.2{t\r) probabilities are 
given in Appendix B. 

Fully stochastic model: An alternative way of account- 
ing for the effects of the full mesoscopic evolution is to 
pre-compute the probability distributions of times for 
completion of the reaction, Pr{t), and binding failure, 
Pf{t). To compute these probabilities, an ensemble of 
trajectories that start at a uniformly chosen position on 
the inner sphere at radius ri is evolved until either the 
substrate binds and reacts or the unbound substrate es- 
capes and passes through an absorbing sphere at inter- 
mediate distance r^ from the binding site. The binding 
probability can be estimated from the fraction of reactive 



trajectories and the probability densities Pr{t) and P/(t) 
can be constructed using analytical fits to the estimated 
cumulative distribution functions obtained from the re- 
action and failure timefp21_ Given this information, once 
the substrate is at ri in step (4), the binding probability 
can be used to determine if reaction will occur and the 
reaction time can be drawn from P^ (t) and used to com- 
plete the cycle, or if no reaction occurs the time can be 
drawn from Pf(t) and used to increment the time. 



IV. SIMULATION OF PGK ENZYME KINETICS 

The simulations employing hybrid MD-MPC dynamics 
were carried out on a system comprising a single PGK 
enzyme with bound ADP, a bPG substrate molecule and 
solvent molecules in a cubic box of length L with peri- 
odic boundary conditions. The units used in the sim- 
ulation are given in terms of length £, mass to, energy 
e and time r. In these units the simulation box had 
length L — AO and contained 640, 000 solvent particles 
of mass m ~ 1, resulting in a density p — 10. The 
mass of the beads comprising the enzyme was taken to 
be M = 10, so that the mass ratio of solvent to beads 
was set to /i = M/m — 10. The solvent particles interact 
with all beads through the truncated repulsive poten- 
tial in Eq. (18 1 with an adjustable a, usually taken to 
be a 



1. Simulation of the enzyme-substrate system 
consists of numerically integrating Newton's equations 
of motion for all bead and solvent particles that inter- 
act with a time step of At — 0.005 for time intervals 
T = 1 between multiparticlc collisions. Information from 
such direct simulations of the dynamics is required for 
both the diffusive encounters between the enzyme and 
substrate and the subsequent binding and reaction pro- 
cesses. These two aspects are discussed in the following 
subsections. 



A. Diffusive dynamics 



long time tails in the velocity correlation function which 
make important contributions to the diffusion coefficient. 
For this reason it is convenient to estimate D by extrap- 
olation of the time-dependent diffusion coefficient to in- 
finite time since 



Dit) 



dt (V(t) • V> - £> 






(6) 



where a^ = {2/'i){A'K{-q + D)y^/'^{mpY/'^ with 77 the 
shear viscosity. The power-law behavior of this quantity 
arises from coupling of the substrate to hydrodynamic 
modes of the solvent. The time-dependent diffusion co- 
efficient is plotted versus i~^/^ in Fig. |4| and shows the 
long-time power-law behavior. For a substrate with mass 
M = 10 in a solvent with p = 10, /c^T = 1/3, substrate- 
solvent Lennard- Jones parameters cr = 0.5 and e = 1, 
we find D = 0.063. Hydrodynamic effects dominate the 
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FIG. 4: The simulated value of the diffusion coefficient com- 
pared to the estimated time-dependent diflusion coefficient, 
D{t), in Eq.([6|, versus t~^''^ for an isolated Brownian par- 
ticle with mass ratio = 10, p = 10, fcsT = 1/3, a = 0.5. 
From the fit of the data, the value of the diffusion coefficient 
is D = 0.063. 



Although the diffusive encounters between the sub- 
strate and enzyme are treated analytically, these calcula- 
tions require the diffusion coefficient D of the substrate 
as input into the analytical formulas. Therefore, in this 
subsection we present results for D for the explicit in- 
teraction and penetrating solvent models. Since the sub- 
strate does not interact with the enzyme in this regime 
we need only consider the motion of the substrate in pure 
solvent. 

Explicit interaction model: In the explicit interaction 
model the substrate interacts with the solvent molecules 
through repulsive Lennard- Jones potentials and the sol- 
vent molecules undergo multiparticlc collisions. The dif- 
fusion coefficient may be determined directly by simu- 
lation from the velocity autocorrelation function or the 
mean square displacement. Hydrodynamic effects are in- 
cluded in the MD-MPC dynamics and these give rise to 



contributions to the diffusion coefficient and it is only 
weakly dependent on the mass of the substrate and sol- 
vent molecules. For a very large substrate molecule the 
diffusion coefficient takes a Stokes-Einstein form and is 
independent of the mass. 

Penetrating solvent model: The diffusion coefficient 
can be computed analytically for the penetrating sol- 
vent model. In the collision step, the rotation matrix 
is uniformly selected from a set of matrices in which the 
rotation by the angles a and —a around a given set of 
axes are equally probable. The operation of the rotation 
matrix on a general vector r for a rotation by angle a 
around a unit vector n can be written succinctly as 

ijj ■ r = r cos a + h{h ■ r){l — cos a) + (r x fi) sin a. (7) 

Since the substrate bead behaves as a point particle with 
respect to hydrodynamic flow, the only contribution to 



the self-difFusion coefficient comes from tfie rotation colli- 
sion step. Hence for this system, the decay of the velocity 
autocorrelation function for an isolated bead is expected 
to be a single exponential. 

The self-difFusion coefficient for this model can be com- 
puted from the velocity autocorrelation function using 
the trapezoidal rule, 



where 



D = 



1 



3 I 2 



dt {V ■ V{t)) 



(8) 



V-V)+Y,{V-V{nT))], (9) 



where r is the collision time and the brackets (• • ■ ) cor- 
respond to an average over the stochastic realizations 
(choice of rotation matrices) and the equilibrium distri- 
bution of the system. If the matrices are chosen uni- 
formly and the rotation angles a and —a are equally 
probable, then the Markovian dynamics for a given cell 
has the limit distribution 



P{n,rn,Vn-R,V) 



^|^n™K)x-n„,(y), (lo) 



where n is the number of solvent particles in the cell 
containing the tagged particle, Vc is the volume of the 
cell (here taken to be unity) and Ii„i{vn) is the normal- 
ized Maxwell-Boltzmann distribution for a system of n- 
particles at temperature T. Using this form, one finds 
that {V ■V)= ^ksT/M, and 
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{V-V{r)) = —J2{V-iv, + u,,-{V-v,))) 

{V-Uj-iV-v,)), (11) 
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i=l 
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where lJ = X]i=i ^i/''^R ^-nd ur is the total number of 
rotation matrices. Inserting the stationary density in 
Eq. (10), and defining the mass ratio /i = Al/m, one 
gets 
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M(l,2 + /i,-p) 
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(12) 



where M(l, 2 + /i, — p) is Kummer's function of the first 
kincP^ and c-y = TraJ/3. If there is no correlation 
between solvent particles occupying the cell containing 
tagged particles following the collision steps, so that 
{V ■ Vinr)} = (1 - j){V ■ V{{n - l)r)), we conclude 
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■pM(l,2 + /i,-p). 



(14) 



The self-diffusion coefficient is plotted in Fig. [5] as a 
function of the mass ratio /i for the simulation values 
ksT = 1/3 and p = 10 and c^ == 1/3. 
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FIG. 5: Plot of the self-diffusion coefficient D for a tagged 
particle in the penetrating solvent model as a function of the 
mass ratio /i. 

While the diffusion coefficient depends weakly on the 
mass ratio for the explicit solvent interaction model, it 
does depends strongly on the mass ratio for the penetrat- 
ing solvent model. In order to facilitate comparisons be- 
tween these two solvent interaction models, we choose the 
mass ratio so that the self-diffusion coefficient of an iso- 
lated bead matches that in the interacting solvent model. 
Note that for the mass ratio /i = 10 used in the inter- 
acting solvent model, the self-diffusion coefficient in the 
penetrating solvent model is substantially larger than in 
the interacting model {D = 0.086 > 0.063), and a mass 
ratio of roughly /i = 28.5 must be used for the dynamics 
of the tagged particle to be comparable. Simulations of 
an isolated Brownian particle immersed in the penetrat- 
ing solvent validate the predictions of Eq. (13 1. Finally, 



we note that the penetrating solvent model without hy- 
drodynamic interactions is also given by Eq. (13). 



B. Substrate binding and reaction 

The position of the bPG substrate was randomly cho- 
sen on a spherical shell at a distance r\ = 1 from the 
active binding site of the enzyme. The distance was cho- 
sen so that the bPG substrate does not interact with the 
active site or other parts of the enzyme. For each real- 
ization of the dynamics, the enzyme configuration was 
equilibrated in the presence of the solvent while con- 
straining the bPG substrate in position. The run was 



then initiated by randomly drawing the bPG velocity 
from a Maxwell-Boltzniann distribution at an effective 
temperature of k^T = 1/3 and releasing the constraint. 
If the substrate bound to the enzyme (determined by a 
distance criterion), the time of binding of the bPG was 
recorded. If instead the distance of the substrate to the 
active site reached a large value, here taken to be at a 
substrate-active site distance of r^ = 9, the evolution of 
a realization was terminated and the failure time was 
recorded. Upon binding, the form of the network po- 
tential for the enzyme allows the enzyme to close to an 
activated form. The time of closing, again determined 
by a distance criterion between conserved, rigid sections 
of the enzyme, was recorded. Once the enzyme closed, 
a reaction time r^ was drawn from a Poisson distribu- 
tion (here taken to have a mean reaction time of Tr = 25 
time units), which defines the rate at which an unbinding 
potential was activated by the control parameter ^. 

The probability densities for the time of substrate 
binding, the closing time of the enzyme after binding, 
and the overall cycle time are shown in Fig. [6J The an- 
alytical fit to the densities with bootstrap estimates for 
uncertainties were computed from the raw data using the 
procedure described in Ref. f44J . A prominent feature in 
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FIG. 6: Probability densities for the full solvent model as a 
function of the collision time. The results are for simulation 
conditions ^ = 10, ksT = 1/3, p = 10, with a solvent-bead 
interaction a = 0.5 cell lengths, corresponding to cr = 2.5A. 

the probability density of binding times is the long alge- 
braic tail, which is a signature of the substrate initially 
moving away from the enzyme but eventually diffusing 
into the active site. The form of the tail in this density is 
consistent with the asymptotic long time behavior for a 
particle diffusing into an absorbing region in three dimen- 
sions. Note that the probability density for the overall 
cycle time can be decomposed into a convolution of the 
density for binding, closing and diffusion away from the 
binding site after the reaction is complete. Since diffu- 
sive motion leads to densities with heavy tails, the overall 



cycle time density is broad, which is characteristic of al- 
gebraic tails. 

Another important qualitative feature of the solvent- 
enzyme model is the variable degree of solvation of the 
bPG substrate during the binding process. When the 
distance a = 0.5 characterizing the solvent-bead repul- 
sion is large enough, the solvent is unable to penetrate 
the volume occupied by the enzyme. The bPG substrate 
binds to a region inside the enzyme that is exposed when 
the enzyme is in an open conformation. Upon binding, 
the enzyme closes via a hinge-like mechanism and brings 
the ADP-bPG substrates near one another enabling the 
transfer of the phosphoryl group. Less solvent is able to 
penetrate into the binding pocket of the bPG substrate 
in the closed conformation of the enzyme, and hence sol- 
vent is expelled from the pocket as bPG binds and the 
enzyme closes , providing a favorable environment for the 
catalysis.liSIlll 

The expulsion of solvent can be tracked by computing 
the local solvent density around the bPG substrate as it 
binds and reacts, as can be seen in Fig. [7J This drying 
effect is highly sensitive to the choice of the repulsive in- 
teraction parameter a. When a = 0.5 (see top panel of 
Fig. I?! , the bound substrate typically has 2 fewer solvent 
particles solvating it, whereas away from the enzyme the 
average number of solvating fluid particles corresponds 
to the value of the bulk density (p = 10). This differ- 
ence between bulk and bound solvation levels increases 
as the repulsion parameter a increases (see bottom panel 
of Fig.lTlwhere a = 0.7). There are important differences 
in the qualitative nature of the dynamics when the repul- 
sion parameter becomes large. Although the exterior of 
the enzyme experiences a larger overall friction, the dis- 
sipating effect of the solvent on the enzyme-substrate in- 
teraction is decreased in the pocket of the enzyme where 
the binding occurs. The bPG substrate retains a high ki- 
netic energy upon entering the pocket for a longer period 
of time due to a limitation in the simple model of the 
binding process in which the substrate effectively inter- 
acts with only a few beads of the enzyme. Because of the 
limited coupling of the beads in the active site to other 
beads in the protein, the excess energy of the substrate is 
slowly dispersed into internal motions of the protein and 
solvent. For this reason, we focus primarily on a regime 
in which the solvent rapidly dissipates energy {a = 0.5). 

In Fig. [8] the probability densities for the binding time, 
enzyme closing time and overall cycle time are presented. 
Looking at the top panel, we see that the probability den- 
sities of the binding time for the interacting and pene- 
trating solvent models are comparable once the dynamics 
has been properly scaled by the mass ratio. This simi- 
larity is not surprising, as the time scale for binding is 
primarily determined by diffusive motion and is not sen- 
sitive to the level of solvation of the substrate by the 
fluid particles. However the absence of hydrodynamic 
flow around the enzyme and substrate has a profound 
effect on both the form of the probability density, which 
is significantly broadened, and the mean binding time. 
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FIG. 8: Probability densities P{t) for substrate binding (top 
panel), enzyme closing (middle panel), and total reaction cy- 
cle (bottom panel) versus time. The black curves correspond 
to results for the interacting solvent model, the red curves 
correspond to the results for the penetrating solvent model 
with hydrodynamics and the blue curves are the results for 
the penetrating solvent model without hydrodynamics. 



namics is significant, and shortens the time required for 
the enzyme to close. 

The overall cycle time density is a convolution of the 
binding time and closing time densities, and is therefore 
different for all three models. 



FIG. 7: Time series showing the reduction in the number of 
solvent particles in the vicinity of the bPG substrate as it 
binds to the enzyme. The red curves show the number of 
solvent particles in the cell containing the bPG substrate as 
a function of time, while the black curves denote the distance 
of the substrate to the enzyme binding site (measured in cell 
length units, where 1 cell length is 5A). (top) a = 0.5, (bot- 
tom) a — 0.7. 



which is shifted by a factor of roughly a factor of three. 
In addition, the binding probability is significantly re- 
duced from Pr — 0.078, in the presence of hydrodynam- 
ics, to Pr = 0.03, which can have a significant impact on 
the density for the overall substrate conversion time when 
the concentration of substrates is elevated. Note that the 
probability density of binding times has a strong tail for 
all models, indicative of the importance of the diffusive 
dynamics experienced by the substrate. 

The time required for the enzyme to close after binding 
is noticeably different in all three models. The penetrat- 
ing solvent model does not account for solvent expulsion 
as the enzyme closes, and therefore has a higher net fric- 
tion and longer time scale than is present in the explicit 
interaction model. Once again, the effect of hydrody- 



C. Fully stochastic model 

A stochastic procedure can be implemented for the 
overall enzymatic process using data from the numeri- 
cal simulations and the computed values of the binding 
probability starting from a radial distance of ri. If the 
binding is accepted starting from the inner sphere with 
probability Pr, which for the explicit solvent model is 
approximately Pr = 0.078, the overall cycle time for the 
reactive process can be added to the overall time for the 
process by drawing from the numerically-obtained prob- 
ability densities and cumulative distributions. To carry 
out the procedure, the reaction time is drawn by numer- 
ically solving the equation Ccycie(iu) = u for the time i„ 
using bisection or Newton-Raphson methods, where u is 
a random variable drawn uniformly from the unit inter- 
val. Here, Ccycic(i) is the cumulative distribution for the 
cycle obtained from the simulation. 

To convert the system collision time into physical 
units, note that the self-diffusion coefficient in system 
units is 0.06 £'^/t. Equating this with the desired value 
of the diffusion coefficient in the cytoplasm of roughly 
D = 1 ■ lO-'^cmVs, we conclude that r = 1.5 • 10"^" 
seconds. Using this scaling, we find that the typical time 
required for the PGK enzyme to close following binding 
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of both substrates is on the order of 3 to 6 ns for the sol- 
vent models incorporating hydrod ynamic flow, which is 
consistent with experimentaPSHSl and simulatiorP^I stud- 
ies of the enzyme domain motions. 

The probability density Pconvit) of substrate conver- 
sion times is shown in Fig. [9] Somewhat surprisingly, 
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FIG. 9: Probability density Pconv(t) of the substrate conver- 
sion time to products versus time expressed in milliseconds for 
the explicit solvent model. The other models yield essentially 
identical results since the substrate conversion is determined 
primarily by diffusion when the substrate is at physiological 
concentrations 

no difference in the probability density of substrate con- 
version time is readily observable at the enzyme concen- 
tration studies here even though the binding probability 
is more than two times larger in the presence of hydro- 
dynamics than in its absence. This is due to the multi- 
ple convolutions of the first passage time densities which 
have heavy and prominent tails that tend to smooth out 
observable differences after multiple convolutions. 



V. SUMMARY 

A stochastic method for computing the probability 
density of the time required for the enzymatic catalysis 
of a substrate to product was constructed. The method 
consists of combining analytical computations of bind- 
ing probabilities and first-passage times of a substrate 
diffusing between two concentric absorbing spheres with 
explicit simulation of motion of the substrate in the im- 
mediate vicinity of the enzyme. Once the explicit sim- 
ulations have been performed and the data analyzed in 
terms of binding probabilities and first passage time den- 
sities, the method allows the probability density of the 
time required for the phosphate transfer to be computed 
at a variety of enzyme concentrations. 

The method was illustrated by considering the cat- 
alytic transfer of a phosphate group from bPG to a 



bound ADP substrate by the phosphoglycerate kinase en- 
zyme under physiological conditions. The binding prob- 
ability and phosphoryl group transfer times for a sub- 
strate diffusing in a 0.1 mM concentration of phospho- 
glycerate kinase were computed under three different 
solvent conditions using a network model of the enzy- 
matic system constructed from the morphing analysis of 
the conformational change between the open and closed 
conformation^28129j qJ ^j-^g enzyme. The solvent models 
were chosen to selectively account for various degrees 
of correlated solvent motion to probe the importance 
of collective flow effects on the enzyme dynamics. It 
was demonstrated that dynamical solvent flow effects as- 
sist the binding of the substrate to the active site of 
the enzyme and facilitate the hinge motion of the en- 
zyme that leads to its closing. Two different models that 
incorporate hydrodynamic flow effects, one with direct 
solute-solvent interactions and another penetrating sol- 
vent model where solvent particles are treated as point 
particles in their interactions with the substrate and pro- 
tein, have similar binding probabilities and cycle time 
densities. However, the density profiles of the solvent 
near the active site as the enzyme closes post-substrate 
binding differ, since expulsion of the solvent from the 
binding pocket is not possible for the penetrating-solvent 
model. In contrast, a Smoluchowski-type model in which 
all beads feel a friction that is independent of the con- 
formation of the enzyme is characterized by a lower sub- 
strate binding probability and a shift in the cycle time 
density to larger time scales relative to the models in- 
corporating hydrodynamic effects. The lower substrate 
binding probability leads to a detectable shift in the max- 
imum appearing in the density of substrate conversion 
times. 

The validity of the stochastic method presented here 
relies on a number of assumptions that are questionable 
for the behavior of the enzymatic system in a cellular 
environment. It has been assumed that the enzymes are 
homogeneously distributed with no correlation between 
their positions in the volume. It is quite possible that the 
enzymes are, in fact, locally clustered in the cytoplasm in 
a way that effectively reduces the distance between them 
and the substrates thereby enhancing their efficiency. 
This is likely to be the case if there is correlation be- 
tween the spatial location of the phosphoglycerate kinase 
enzyme and enzymes such as glyceraldehyde phosphate 
dehydrogenase that act earlier in glycolysis. In addition, 
it has been assumed that the dynamics of the substrate in 
the complex, crowded cytoplasm is diffusive, which may 
be reasonable on long time scales but less accurate on the 
time scale of solvent motion. However, subdiffusive mo- 
tion of proteins and finite-size probe molecules has been 
seen in crowded cellular environments."*^ —Nonetheless, 
assuming substrates do move diffusively in the cytoplasm 
at long times, the diffusive nature of the substrate dy- 
namics leads to a broad distribution of substrate conver- 
sion times that differs substantially for the exponential 
distribution one might anticipate from mass action kinet- 



11 



It is straightforward, though coraputationaUy inten- 
sive, to incorporate more detailed models of the en- 
zymatic system to produce quantitatively accurate re- 
sults. This is readily accomplished by performing all 
atom simulations of the system complete with detailed 
molecular mechanical-based interaction potentials and 
quantum-mechanical analysis of chemical reaction path- 
ways. Nonetheless, it is likely that the observation that 
the solvent flow assists the binding and subsequent pro- 
tein motions will also be observed in more detailed mod- 
els of the enzymatic system. 
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Appendix A: PGK potential functions 

In this Appendix we give the detailed form of the po- 
tential function Vps that governs the dynamics of the 
protein and its interactions with the bPG substrate. 

Bonds in the set Be of common links were assigned 
bond potentials Vc{rij) constructed in the following way. 
The potentials for the common links in the open and 
closed configurations of the enzyme, Vco and Vcc, are 
given by 



Vcc 



kh 



E ( 






<ij>eBha 



(15) 



E 



<ij>eBsc 




where the parameters l]""^' and a^ 



i°''^ and (t;°'^ were determined 
by the equilibrium distances for the harmonic and soft- 
common links in the open and closed conformations and 
k}i is the force constant for the hard elastic network 
bonds. Given this input, the potential for the common 
interactions Vc was taken to be the lowest eigenvalue of 
a two-dimensional empirical valence bond (EVB) matrix 
with constant off-diagonal elements A, so that^^ 

K = ^ {{Vco + Vcc) - {{Vco - Vcc)^ + 4A2) '/') . (16) 

This form of the potential allows the system to smoothly 
switch between stable open and closed configurations. 



Links in the soft-open, Bso, and soft-closed, Bsx sets were 
assigned bond potentials 



VsiR^J)^e{5 
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(17) 



with identical forms. In addition, monomeric beads rep- 
resenting amino acid residues repel one another at short 
distances according to a truncated Lennard- Jones (LJ) 
potential 



Vr 




<ij> 



Kij 



+ 1 0{abb ~ Rij), 

(18) 

where 6{x) is the Heaviside function. The bPG substrate, 
represented by a single bead with coordinate R, also in- 
teracts with all beads in the protein through a repulsive 
LJ potential of this form, Vy (i?&i), where Rbi = |R— R^j 
with ebs and at,s energy and distance parameters. 



Interactions governing the reactive event and 
conformational changes 

The binding interaction F/^^(R,Rg, RJ,R5) between 
the bPG substrate at position R and the enzyme was 
designed to depend on the distance between the sub- 
strate and bead with coordinate R", as well as the ori- 
entation of the substrate with respect to a coordinate 
frame determined by three beads defining the binding 
pocket of the enzyme. Defining the relative position vec- 
tor Rsi — R — RI = RsiR with magnitude Rsi and 
direction R of the substrate with respect to a coordinate 
system centered on the binding site i?°, the projection 
Rgi — R ■ {RIq X ^21 ) is computed, where Rf^ is the 
unit vector along i?° = Rf — R'-. The binding potential 
is then taken to be 
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(19) 



where Ks = 1.5 in the energy units. In Eq. (19), J{R) is 
a smooth cut-off function 



I{R) 



1, 

(R^~R? (r, _ on 

0, 



R<Ri 

2R) , Ri ^ R ^ Ru 

R > Ru 



(20) 
where the upper and lower cut-off values are set to 
Ru — 3cr and Ri — 2.5cr. The potential insures that the 
optimal angle of approach and binding of the substrate 
in the active site pocket is along the J?2i ^ -^lo direc- 
tion. In principle, the excluded volume interactions of 
the substrate bead with the enzyme beads are sufficient 
to determine the binding pathway of the substrate, while 



12 



the orientational dependence of the binding potential in 
Eq. ( 19 ) restricts the binding location in the active site. 
As the substrate binds it triggers conformational 
changes in the protein that lead to hinge closing to bring 
the bPG and ADP substrates into proximity for the phos- 
phoryl group transfer. Thus, as bPG interacts with the 
protein in the course of binding to the active site, the 
open protein configuration is destabilized with respect 
to the closed configurations, driving the enzyme towards 
the closed conformation. To achieve this conformational 
change in the network model, the interaction potentials 
for the soft, non-common set of links are modified. We 
define the reaction coordinate ^, where 



^ = - (1 + tanhx) , 



where 






fci 



Rii)' 



{Rbi - i?gj^ 



(21) 



(22) 



and R^^ is the initially large distance between the sub- 
strate and the binding site in the enzyme in the open 
configuration and R^^ is the same distance in the bound, 
closed complex. Since the substrate is unbound and 
hence far from the binding pocket in the starting con- 
figuration, R'^j^ ^ i?gj^. Note that when the substrate is 
far from the enzyme, x is large and negative and ^ ~ 0, 
whereas x becomes large and positive as the substrate 
moves towards the binding site with the result that ^ « 1 
upon binding. Given this reaction coordinate, the soft, 
non-common potential function is taken to be 

<ij>eBs^ <ij>eBsa 

The protein-substrate interaction potential is given by 
the sum of these contributions: 

Fps(R'^^R;^(i^5l)) = V^c+l^r+K^'^'+V^'^'+Kc (24) 

After binding, the reaction coordinate ^ is treated as 
an external control parameter that is governed by the 
equation: 



m 



t/Xr li t < Tr 

otherwise, 



(25) 



where r is the reaction time drawn from an exponential 
distribution P{tj.) = r^^e"'^''''^'' and r^ is the average 
reaction time. Upon completion of the reaction when 
^ = 0, the interaction between the substrate in the bind- 
ing pocket and the binding site is changed to a repulsive 
Lcnnard- Jones interaction to reflect the unstable inter- 
action of the altered substrate and the binding pocket. 
Since ^ = 0, the closed configuration is unstable and the 
enzyme reopens, completing the cycle. In this treatment, 
the reaction is treated irreversibly and the surrounding 
solvent absorbs energy from the chemical process, leading 



to a slight heating of the solvent. The average reaction 
time Tr is taken to be 25 time units, corresponding to a 
physical reaction time of roughly 2.5 ns. Note that the 
precise value for the average reaction time is unimpor- 
tant for looking at the qualitative effects of the solvent 
environment on the dynamics of the enzymatic system. 
Detailed quantum chemical calculations are required to 
determine if this estimate of the reaction time from the 
metastable bound state to a final state consisting of the 
products bound in a closed conformation of the enzyme 
is reasonable. 



Appendix B: Diffusion of substrate to a region near 
enzyme 

The first passage time distribution P(f |ro) for a Brow- 
nian walker starting from position Tq at time i = onto 
a sphere centered at the origin can be computed from the 
survival probability distribution F(i|ro) using 



P{Ar^) = - 



dF{t\ro) 



dt 



(26) 



sphere 



where the derivative of F{t) only includes the flux of 
walkers into the sphere and 



P(i|ro)= f drPir,t;ro). 



(27) 



In Eq. (27 1, P(r,i;ro) is the conditional probability of 



finding the walker at position r at time t given that it 
was initially at rg, and Q is the domain of the system. 
We shall assume that the walker is confined between two 
absorbing spheres of radii r_ = ri and r+ — r2. Given 
the spherical boundaries of the domain, it is natural to 
express positions in terms of spherical polar coordinates 
(r, 9, (p), where the z-axis from which the angle 9 is mea- 
sured relative to the vector connecting the origin to a 
specific point r^ on the inner sphere. The angle (p can 
measured from the plane containing the vectors Tp and 
Tq so that 00 = 0. The evolution of the conditional prob- 
ability is determined by the diffusion equation 



dP 



DViP, 



(28) 



and satisfies the boundary condition P(r,0;ro) = 6{r — 
rg). From the diffusion equation, we find that the first- 
passage distribution through a spherical domain at radial 
distance r_ is given by 



= -D / drViP 
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(29) 
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where the second Hne follows from Green's theorem. The 
domain il contains all points with radii in the range 
[r_,r_|_], and the integral over the inner sphere can be 
written to obtain 



P{t\ro) = Dr^_ 



ddi 



dP{r^,9,cl),t;ro) 
dr 



(30) 

The diffusion equation may be solved for arbitrary co- 
ordinates r and ro in the presence of absorbing bound- 
aries by expanding the density P(r, i|ro) in spherical po- 
lar coordinates. The absorbing boundary conditions re- 
quire that 



P(r_, 0,0,t;ro, 00, 0o)=O 
P(r+, 0,0,t;ro, 00, <^o)=0. 



(31) 




Although a general series solution in spherical har- 
monic functions for P(r,t|ro) is possible, the spheri- 
cally averaged flux F[t\ro) and first-passage time dis- 
tribution P(i|ro) are simple to obtain since only the 
first, spherically-symmetric term in the expansion re- 
mains. From the differential equation for the expan- 
sion coefficients, one finds that the Laplace transform 
P(s|ro) — J^ dte^^*P{t\ro) of the first passage time 
density P(t|ro) for the inner sphere is given by^ 



P{s\ro) 



1/2 



Cl/2(a;o,a;^ 
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xo J Ci/2ix-,x+) 
r\ ^1^ Ci/2(a;o,a:+) 
tq) Cii2{x-,x+y 



(32) 
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where xq is the scaled variable xo = \J s/ Dr^), x^ = 
\/s7^r+, X- = yJs/Dr_, and Cu{a,h) = I^{a)K^{b) - 
Iiy{b)K^{a), where Iv{x) and K^{x) are modified Bessel 
functions. For a large outer sphere for which r_|_ ^ r_, 
Ci/2{x,x+) — >■ —Ii/2(x+)Ki/2{x). Considering a particle 
that can start at any point on a spherical shell at r = ro, 
we can write P(s|ro) -- ^/rZJrQ Ki/2{xq)/Ki/2{x^). 
Noting that 



FIG. 10; Absorption time probability density versus time. 
The top panel is the absorption time for the absorption onto 
an inner sphere at ri = 7 starting from a radial distance r = 
10 in length units l. The bottom panel shows the absorption 
time density (top) and cumulative distribution (bottom) for 
the outer sphere, where the outer absorbing sphere radius is 
set to be r2 = 31.6 and r = 10. 



fco(x) 



"" Ki/2{x) = ^e 



the Laplace transform Pi{s\ro) of the first passage den- 
sity to the inner sphere can be approximated by 



A(s|ro) 



kojxo) 
fco(x_) 



/s/D (ro-r_) 



(33) 



which can be explicitly inverted to obtain the normal- 
ized first-passage distribution Pi{t\r) for particles that 
are absorbed at the inner sphere radial distance ri start- 
ing from the spherical shell at distance r, 



Pi(t|r) 



V4:7rDt^ ' 



-(r-riy/{4Dt) 



(34) 



This result is plotted in Fig. 10 (top panel). Note that 



the fraction of particles absorbed at the inner sphere in 
the infinite time limit can be computed from the s = 
limit of Eq. (32 1, yielding 



Pi(s = 0|r) = Pi(r) = 



n r2 



r r2 - ri 



The fraction of particles absorbing at the outer bound- 
ary in the infinite time limit is P2(?') — ^^ Pi{f)- These 
probabilities play an important role in the stochastic sim- 
ulation algorithm. 

The first-passage time density at the outer sphere is 
obtained similarly, although the inversion of the Laplace 



transform P2(s|r) is complicated since 



P2is\r) 



r) Ci/2{xi,X2) 
r2 sinhxi e^^ — sinha; e~^i 
r sinhxi e^^^ — sinhx2 e^^i 



where x = ^J sjDr and Xi = ^Js/Dri. Although the 
density can be approximated using series expansions for 
0-functions, it is a simple matter to inver t P2{s\r) nu- 
merically using the Stehfest algorithrrPlMl, 

To draw a random time from the first-passage den- 
sity Pi(i|r), one first defines the cumulative distribution 
Ci{t\r) = /odrFi(T|r) = 1 - erf((r - ri)/\/4Di). Sup- 
pose u is drawn uniformly from the unit interval. Setting 
u = Ci{tu\r) and solving for t„ gives 



{r -rif 



AD (icrf(l-M))' 



(35) 
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where ierf is the inverse error function which can be 
solved for numerically in an efficient manner using the 
secant method. The set i„ are then drawn from the first- 
passage distribution. 



The task of drawing from the distribution P2{t\r) 



shown in the bottom panel of Fig. 10 is readily accom- 
plished by drawing a random number p uniformly on 
(0, 1) and then solving the implicit equation C2(iu|r) = p 
for the time i„, where C2{t\r) is the cumulative distribu- 
tion C2{t\r) — ^^dr P2{T\r). The cumulative distribu- 
tion can be computed numerically by applying the Ste- 
hfest algorithm to form the inverse Laplace transfo rm o f 
C2{s\r) = P2{s\r)/s (see bottom-most panel of Fig. 10). 
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